#include "nbodies_final.h"

//TODO faire un script  pour pouvoir le lancer (qui appelle le script générant le fichier)

// calcule la distance entre deux masses
float dst(struct body * a, float x, float y)
{
  return sqrt(((x-a->x)*(x-a->x) + (y-a->y)*(y-a->y))); 
}

//calcule une des composantes du vecteur unitaire entre deux masses
float unit(int dir, struct body * a, float x, float y)
{
  if(dir == X)
    {
      float length = (x - a->x);
      return (length/dst(a,x,y));
    }
  else
    {
      float length = (y - a->y);
      return (length/dst(a,x,y));
    }
}

//Fonction de debug
void print(struct body ** bodies, int nb_bodies)
{
  int i;
  for(i=0 ; i<nb_bodies ; i++)
    {
      printf("%d %f %f\n",i,bodies[i]->x,bodies[i]->y);
    }
  printf("\n");
}

//Fonction de récupération des données initiales en fonction du rang.
struct body ** parse(char * path,int * nb_bodies,int rank)
{
  FILE * fd = fopen(path,"r");
  int nb_all;
  fscanf(fd,"%d %d",&nb_all,nb_bodies);
  struct body ** bodies = malloc((*nb_bodies)*sizeof(struct body *));
  int i=0,j;
  float* ligne = malloc(5*sizeof(float));
  while(i<(rank*(*nb_bodies)))
    {
      fscanf(fd,"%f %f %f %f %f",&ligne[0],&ligne[1],&ligne[2],&ligne[3],&ligne[4]);
      i++;
    }
  
  for (j=0 ; j<*nb_bodies ; j++){
    bodies[j] = malloc(sizeof(struct body));
    bodies[j]->speed = malloc(sizeof(struct speed));
    bodies[j]->new_acc = malloc(sizeof(struct acc));
    bodies[j]->new_acc->x = 0.0;
    bodies[j]->new_acc->y = 0.0;
    fscanf(fd,"%f %f %f %f %f\n",&bodies[j]->m,&bodies[j]->x,&bodies[j]->y,
	   &bodies[j]->speed->x,&bodies[j]->speed->y);  
  }
  return bodies;
  fclose(fd);
}

//Fonction de calcul des forces sur les masses
int calc(float *buffer, struct body **bodies, int nb_bodies, int first, int iter)
{
  static int call;
  call+=1;
  int myrank;
  MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
     
  float G = 6.67384*pow(10,-11);
  int i,j;
  if(first==0)//calcul de l'influence des masses du processus appelant
    {
      for(i=0 ; i<nb_bodies ; i++)
	{
	  //calcul des forces sur 1 masse.
	  
	  bodies[i]->new_acc->x = 0;
	  bodies[i]->new_acc->y = 0;
	  bodies[i]->close_m=-1;
	  for (j=0 ; j<nb_bodies ; j++)
	    {
	      if(i!=j)
		{
		  float dist = dst(bodies[i],bodies[j]->x, bodies[j]->y);
		
		  float f = G*bodies[j]->m/(dist*dist);
		  bodies[i]->new_acc->x += f*unit(X,bodies[i],bodies[j]->x,bodies[j]->y);
		  bodies[i]->new_acc->y += f*unit(Y,bodies[i],bodies[j]->x,bodies[j]->y);
		  //maj de la planete la plus proche		  
		  if(bodies[i]->close_m == -1 || (dist < bodies[i]->close_dst)){
		    bodies[i]->close_m=bodies[j]->m;
		    bodies[i]->close_x=bodies[j]->x;
		    bodies[i]->close_y=bodies[j]->y;
		    bodies[i]->close_dst=dist;
		  }
		  		  
		}
	      
	    }
	}
    }
  else //calcul de l'influence des masses des autres processus
    {
      for(i=0 ; i<nb_bodies ; i++)
	{
	  //calcul des forces sur 1 masse.
	  for (j=0 ; j<nb_bodies ; j++)
	    {
	      float dist = dst(bodies[i],
			       buffer[j*3+iter*nb_bodies*3+1], 
			       buffer[j*3+iter*nb_bodies*3+2]);
	      float f = G*(bodies[i]->m*buffer[j*3+iter*nb_bodies*3])/(dist*dist);
	      if(dist == 0){
		printf("i : %d, j : %d, f : %f, dist : %f call : %d\n",
		       i,j,f,dist,call);
		
	      }
	      bodies[i]->new_acc->x += f*unit(X,bodies[i],buffer[j*3+iter*nb_bodies*3+1],
					      buffer[j*3+iter*nb_bodies*3+2]);
	      bodies[i]->new_acc->y += f*unit(Y,bodies[i],buffer[j*3+iter*nb_bodies*3+1],
					      buffer[j*3+iter*nb_bodies*3+2]);
	      //printf("new_acc : x=%f y=%f\n",bodies[i]->new_acc->x,bodies[i]->new_acc->y);
	      //maj de la planete la plus proche
	      if(dist < bodies[i]->close_dst || bodies[i]->close_m == -1){
		bodies[i]->close_m=buffer[j*3+iter*nb_bodies*3];
		bodies[i]->close_x=buffer[j*3+iter*nb_bodies*3+1];
		bodies[i]->close_y=buffer[j*3+iter*nb_bodies*3+2];
		bodies[i]->close_dst=dist;
	      }
	    }
	
	}
    }
  return 0;
}



float maj_dt(struct body ** bodies, int nb_bodies, float old_dt)
{
 
  float delta_t;
  int j;
  float delta;
  float dt;
  float a;
  float b;
  float c;
  //calcul de dt
  for(j=0 ; j<nb_bodies ; j++)
    {
      if(bodies[j]->new_acc->x == 0 && bodies[j]->new_acc->y == 0)
	{
	  dt = old_dt;
	}else{
	a = 0.5*sqrt(bodies[j]->new_acc->x*bodies[j]->new_acc->x+
		     bodies[j]->new_acc->y*bodies[j]->new_acc->y);
	
	b = sqrt(bodies[j]->speed->x*bodies[j]->speed->x+
		 bodies[j]->speed->y*bodies[j]->speed->y);
	
	c = -0.1*dst(bodies[j], bodies[j]->close_x,bodies[j]->close_y);
	
	delta = b*b-4*a*c;
	// on choisit la plus grande racine pour ne pas réduire dt inutilement
	dt = (-b+sqrt(delta))/(2*a);
      }
      if(j==0 || dt < delta_t ) delta_t = dt; 
    }
  return delta_t;
}

//Fonction de mise à jour des positions des masses 
//en fonctions des accélérations appliquées.
int maj(struct body ** bodies, int nb_bodies, float dt) 
{
  int j;
  for(j=0 ; j<nb_bodies ; j++)
    {
      bodies[j]->x += dt*bodies[j]->speed->x + bodies[j]->new_acc->x*dt*dt/2; 
      bodies[j]->y += dt*bodies[j]->speed->y + bodies[j]->new_acc->y*dt*dt/2; 
      bodies[j]->speed->x += (bodies[j]->new_acc->x)*dt;
      bodies[j]->speed->y += (bodies[j]->new_acc->y)*dt;
    }
  return 0;
}

//print dans un fichier des informations:
//flag = 0, on ecrit les positions dans le fichier rank.res
//flag = 1, on ecrit les accélérations dans rank.acc
void printToFile(struct body ** bodies, int nb_bodies, int proc, int flag)
{
  char * path;
  char a[15];
  sprintf(a,"%i",proc);
  if(flag ==0)
    {
      path = strcat(a,".res");
    }  
  else
    {
      path = strcat(a,".acc");
    }
  
  FILE * fd = fopen(path,"w+");
  int i;
  if(flag ==0)
    {
    for(i=0 ; i<nb_bodies ; i++)
      {
	fprintf(fd,"%f %f\n",bodies[i]->x,bodies[i]->y);
      }
  }
  else
    {
      for(i=0 ; i<nb_bodies ; i++)
	{
	  fprintf(fd,"%f %f\n",bodies[i]->new_acc->x,bodies[i]->new_acc->y);
	}
    }
  fprintf(fd,"\n");
  fclose(fd);
}

int main(int argc,char * argv[])
{
    if(argc<4)
    {
      printf("utilisation : %s fichier_entrée nbre_de_tours dt_initial\n",argv[0]);
      return 0;
    }
    // chaque proc a le même nombre de corps (qui divise le nombre total de corps).
    // cf fichier d'entrée avec des masses nulles
    char * path = argv[1];
    int max_run = atoi(argv[2]);

    float old_dt = atof(argv[3]);
    float dt;
    int nb_proc;
    int myrank;
    int tag=42;
    MPI_Status status;
    MPI_Init(NULL,NULL);
    MPI_Comm_size(MPI_COMM_WORLD,&nb_proc);
    MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
    int i,j,k;
    int nb_bodies;
    struct body ** bodies = parse(path,&nb_bodies,myrank);
    //structure du buffer float : m px py
    float * buffer = malloc(2*3*nb_bodies*sizeof(float));
    //initialisation du buffer de floats
    for(j=0 ; j<nb_bodies ; j++)
      {
	buffer[j*3] = bodies[j]->m;
	buffer[j*3+1] = bodies[j]->x;
	buffer[j*3+2] = bodies[j]->y;
      }
    
    MPI_Barrier(MPI_COMM_WORLD);
    for(k=0;k<max_run;k++)
      {   
	int iter = 0;
	//initialisation des comm persistantes.
	MPI_Request odd_request,even_request;
	MPI_Request recv_request;
	int size;
	MPI_Pack_size(3*nb_bodies,MPI_FLOAT,MPI_COMM_WORLD,&size);
	int Bbuff[2*size+2*MPI_BSEND_OVERHEAD];
	MPI_Buffer_attach(Bbuff,2*size+2*MPI_BSEND_OVERHEAD);
	MPI_Bsend_init(buffer,3*nb_bodies,MPI_FLOAT,(myrank+1)%nb_proc,tag, MPI_COMM_WORLD,&even_request);
	MPI_Bsend_init(&buffer[3*nb_bodies],3*nb_bodies,MPI_FLOAT,(myrank+1)%nb_proc,tag, MPI_COMM_WORLD,&odd_request);

	for(i=0;i<nb_proc-1;i++)
	  {
	    //--------Communications bloquantes non persistantes
	/*     MPI_Sendrecv(&buffer[iter*3*nb_bodies], 3*nb_bodies, MPI_FLOAT, (myrank+1)%nb_proc, tag, */
/* 			 &buffer[(1-iter)*3*nb_bodies], 3*nb_bodies, MPI_FLOAT, (myrank-1+nb_proc)%nb_proc, tag, */
/* 			 MPI_COMM_WORLD, &status); */
	    //-------------------
	    //------Communications persistantes:
	    if(iter==0)
	      {
		MPI_Start(&even_request);
	      }
	    else
	      {
		MPI_Start(&odd_request);
	      }
	    //------Receive
	    MPI_Irecv(&buffer[(1-iter)*3*nb_bodies],3*nb_bodies,MPI_FLOAT,(myrank-1+nb_proc)%nb_proc,tag,MPI_COMM_WORLD,&recv_request);
	    //------Calcul

	    //-----------------------
	    calc(buffer, bodies, nb_bodies, i, iter);
	    
	    if(iter==0)
	      {
		MPI_Wait(&even_request,&status);
	      }
	    else
	      {
		MPI_Wait(&odd_request,&status);
	      }
	    MPI_Wait(&recv_request,&status);
	    iter = 1 - iter;
	  }
	//dernier calcul
	calc(buffer, bodies, nb_bodies, i, iter);
	
	//maj du dt optimal local
	dt = maj_dt(bodies, nb_bodies, old_dt);
	
	//Libération des requetes
	MPI_Request_free(&odd_request);
	MPI_Request_free(&even_request);

	//mise en commun des dt
	//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	float new_dt;
	MPI_Allreduce(&dt, &new_dt, 1, MPI_FLOAT, MPI_MIN,
				   MPI_COMM_WORLD);
	//	printf("reduce : %d\n",reduce);
	dt = new_dt;
	if(myrank == 0) printf("dt : %f\n",dt);
	//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	
	//mise à jour des positions et vitesses 
	//en fonction des accélérations stockées et du dt commun.	
	maj(bodies, nb_bodies, dt); 
	
	//Print dans un fichier, chaque proc a son fichier.
	printToFile(bodies,nb_bodies,myrank,0);

	// maj buffer 0
	for(j=0 ; j<nb_bodies ; j++)
	  {
	    buffer[j*3] = bodies[j]->m;
	    buffer[j*3+1] = bodies[j]->x;
	    buffer[j*3+2] = bodies[j]->y;
	  }
	void * b;
	int s;
	MPI_Buffer_detach(&b,&s);
      }
    
    for(i=0;i<nb_bodies;i++)
      {
	free(bodies[i]->speed);
	free(bodies[i]->new_acc);
	free(bodies[i]);
      }
    
    free(bodies);
    
    MPI_Finalize();
    return EXIT_SUCCESS;
}
